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Abstract. Estimating the relative alignment between the frontier molecular 
orbitals that dominates the charge transport through single-molecule junctions 
represents a challenge for theory. This requires approaches beyond the widely 
employed framework provided by the density functional theory, wherein the Kohn- 
Sham “orbitals” are treated as if they were real molecular orbitals, which is 
not the case. In this paper, we report results obtained by means of quantum 
chemical calculations, including the EOM-CCSD (equation-of-motion coupled- 
cluster singles and doubles), which is the state-of-the-art of quantum chemistry 
for medium-size molecules like those considered here. These theoretical results 
are validated against data on the molecular orbital energy offset relative to the 
electrodes’ Fermi energy extracted from experiments for junctions based on 4,4’- 
bipyridine and 1,4-dicyanobenzene. 
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1. Introduction 

The last decades marked significant advances in the fabrication and characterization 
of a variety of molecular electronic devices based on a single or a small number 
of molecules. For understanding an impressive amount of experimental material 
accumulated, results of numerical calculations are usually presented. They are 
performed within theoretical frameworks, which are often completely opaque and 
preclude a straightforward interpretation in terms of properties having a simple and 
clear physical meaning. 

An important property of this kind is the relative alignment of the frontier 
orbitals relative to the electrodes’ Fermi energy Ep of a molecule embedded in a 
molecular junction m- The energy offset Eq = min {Epu mo - Ep.Ep - Ehomo) 
of the highest occupied or lowest unoccupied molecular orbital (HOMO or LUMO, 
respectively), whichever is closest to the Fermi level Af, is usually compared to a 
tunneling energy barrier [2]. It is a key quantity in molecular transport, because it 
controls the charge transfer efficiency. Current experimental methods to estimate the 
energy offset Eq employ ultraviolet photoelectron spectroscopy (UPS) [3], thermopower 
Hiniiiiz], and transition voltage spectroscopy (TVS) [8]. TVS is a method that 
became very popular in the molecular electronic community due to its simplicity 
[HI m Uni HH HH HSl im HSl ini HZl- TVS-based results deduced from ref. [5], a 
work that significantly contributed to the TVS popularity, represent an essential piece 
of experimental data to be used in the present study. 

Postulating the existence of a single orbital that dominates the charge transport in 
a molecular junction might appear to be a too crude approximation. However, recent 
extensive analysis [m [m EDI EB Ea EH El] of existing transport data measured 
for a variety of molecular junctions demonstrated that, in the entire voltage accessed 
experimentally, current-voltage {I—V) curves can indeed be excellently be reproduced 
by assuming just the contribution of a single dominant molecular orbital (MO). 
Although this is an enormous simplification, the quantitative description of the relative 
alignment Eq of the dominant orbital of the embedded molecule remains an important 
challenge for ab initio approaches to the charge transport. The vast majority of the 
theoretical approaches of the charge transport through molecular devices utilized to 
date are based on the combination of the nonequilibrium Green’s functions (NEGF) 
and density functional theory (DFT) [D^. The most important drawback of such 
approaches directly related to the main issue considered in this paper is the fact that 
they treat the eigenvalues of the Kohn-Sham equations as if they were energies with 
physical meaning. In reality, as is well known, Kohn-Sham orbitals are mathematical 
objects rather than true molecular orbitals ESj- So, it is not at all surprising that 
their “energies” cannot provide an adequate physical description [DD] EZ] ■ Unoccupied 
orbitals are especially difficult to describe theoretically [221E1]- Therefore, molecular 
junctions exhibiting an n-type (LUMO-mediated) conduction, like the ones to be 
considered in this study, deserve a special consideration. 

Demonstrating that valuable information on the LUMO alignment in molecular 
junctions can be obtained within reliable quantum chemical approaches beyond DFT 
represents an important aim of the present paper. To validate the theoretical approach 
developed here, we will employ experimental data for molecular junctions based on 
1,4-dicyanobenzene (BDGN) [2] and 4,4’-bipyridine (44BPY) [6|. 
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2. Disentangling the contributions to the dominant MO energy offset 

As a remedy of the main drawback of approaches to the charge transport through 
molecular devices that combine nonequilibrium Green’s functions (NEGF) and DFT 
mentioned in Introduction, in more elaborate (so-called “DFT-I-S”) developments 
[MIE] the occupied (cp < Ep) and unoccupied (cp > Ep) KS “orbitals” are rigidly 
shifted in opposite directions by the same (p-independent) amount A obtained either 
by fitting experimental data |29| or in a two-step procedure as follows. First, a value 
^0 = EpjjMO — Ehomo of the HOMO-LUMO gap is deduced from the energies 
E of the various (neutral, anionic and cationic) charge species in the gas phase 
Ao = £anion+Ecation-‘^£neutraU & method known as A-SCF [55] (or, more appropriate, 
A-DFT |30l [311 |22l [27]). Because the value Aq is usually much too large, it is then 
renormalized (Aq —>■ A) by considering image charges of a LUMO (HOMO) modeled 
as a point-like electron (hole) formed in electrodes taken as infinite plates in the 
immediate vicinity of the active molecule [slIEHlE]- Although the renormalization 
found in this way may render the corresponding A-value compatible to experiments, 
recent work |33llMl ESlIlT] has drawn attention on the fact that these assumptions 
are inadequate for realistic molecular junctions |33l|3S[35l[27]. 

An aspect on which we want to draw attention in this study is the following. 
NFGF-DFT transport calculations done as described above utilize an extended 
molecule, which includes several atomic layers from electrodes in addition to the 
active molecule. On the other side, as it has been long recognized, the classical 
image potential originates from the interaction of the electron in the LUMO with 
electronic collective (long-wavelength polarization) modes in the metals, in particular, 
surface, interface and bulk plasmons |36l EZ] EHl |39l HO] [41]. So, the effect of the 
electrons of the atoms of the electrodes belonging to the extended molecule is also 
accounted for in the interaction energy with the image charges. Therefore, the 
procedure of applying E-corrections (these corrections are due to image charges) on 
top of NFGF-DFT approaches is plagued by double counting. 

To overcome this drawback, we propose here a disentangling procedure, which we 
then validate by comparing molecular junctions based on two molecules, namely 1,4- 
dicyanobenzene (BDCN) and 4,4’-bipyridine (44BPY). Previous studies demonstrated 
an n-type (LUMO-mediated) conduction for both types of junctions [U |42l |6l|22] . So, 
it is the LUMO on which attention will be focused below. 

The disentangling scheme for the LUMO energy is presented in Figured) (For the 
numerical values given in Figured] please refer to Sec. |4|) In the absence of molecule- 
electrode couplings, the LUMO energy Ei is given by the lowest electron attachment 
energy EA of the isolated molecule taken with reversed sign {Ei = —EA, Koopman’s 
theorem). The LUMO energy EpuMO of the molecule embedded in a nanojunction 
differs from that of the isolated molecule because an electron transferred to the LUMO 
interacts with the electrodes. The scheme proposed here holds in cases (and we will 
show below that such cases do exist) where it is possible to split this interaction into 
a short-range and a long-range part that can be analyzed separately. They yield two 
contributions (denoted A and ^im, respectively) to the corresponding LUMO energy 
shift 

Elumo = El + A.+ ^ im • (1) 

Eq. (dj) allows one to express the LUMO energy offset as 

£0 = EpuMO — Ep = —EA -I- A -|- — Ep. 


( 2 ) 
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Zero energy (vacuum level) 


<P,„^=-1.37q\J 


E,,,=-EA, =-O.OSeW 


E. .= .EA,= -0.72eV 


0„„2 = -O.69eV 


fo., = r.53eV 


Electrode Fermi energy 


1M4BPY 


2=BDCN 


Figure 1. Disentangling the contributions to the LUMO energy offset for the 
two molecules — 4,4’-bipyridine (44BPY) and 1,4-dicyanobenzene (BDCN) — 
embedded in the molecular junctions considered in the present paper. The electron 
attachment energy with reversed sign represents the LUMO energy of the isolated 
molecule {Ei = —EA). ^im and A represent LUMO energy shifts due to image 
charges and local (contact) effects resulting from long-range and short-range 
interactions between the embedded molecule and electrodes, respectively. The 
nearly equal quantities Ai Pr; A 2 reflect the similar chemical linkage (nitrogen- 
gold affinity) to electrodes. See the main text for details. 


To avoid the double counting issue mentioned above, in addition to the electrodes’ 
collective effect embodied in the image charge contribution we will consider the 
local LUMO energy shift A due to the interaction of the active molecule with the 
metal atoms at the its ends. This accounts for the well known fact that a chemisorbed 
molecule strongly coupled to the substrate often has valence molecular orbital energies 
substantially different from that in the gas phase [43] . Experimental data on molecular 
junctions, indicating a substantial MO energy shift due to local (interface) dipoles |3], 
may also be taken as a confirmation of this hypothesis. 

3. Method 

To obtain the theoretical results presented in this paper we have utilized the 
EOM-CCSD (equation-of-motion coupled-cluster singles and doubles) jH] US] HSj. 
This method represents the state-of-the-art of quantum chemistry for medium-size 
molecules, like the ones to be considered here. The CCSD calculations were performed 
with the CFOUR package [47]. For comparison purposes, results based on hybrid 
coupled clusters (CC2) [48] and regular (strict) second-order algebraic-diagrammatic 
constructions [ADC(2)] |49l [50] will also be presented. ADC(2) calculations have 
been done with the fully parallelized PRICD-I](2) code [5T], which is interfaced to 
MOLCAS [52|. 

The molecular geometries were optimized at the DFT level using the B3LYP 
hybrid functional as implemented in GAUSSIAN 09 [53], a package also employed 
to estimate the electroafhnities by means of A-DFT calculations [26l EQ] [27]. All 
results of the quantum chemical calculations reported here were obtained by employing 
aug-cc-pVDZ (Dunning augmented correlation consistent double zeta) basis sets. As 
shown in recent studies [an m [Ml [27], these basis sets include sufficient diffuse basis 
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functions to properly describe anionic states, and the corresponding results for electron 
attachment energies can be trusted. 

4. Results and discussion 

Obtaining estimates for A from quantum chemical calculations by comparing the 
LUMO energies of an isolated molecule and the same molecule with one or a few 
metal atoms attached at each of its ends will be the object of a further investigation. 
Here we will confine ourselves to quantify the difference in the LUMO energy offsets 
£ 0,2 — £ 0,1 for two molecular junctions consisting of molecules (labeled 1 and 2) joined 
to electrodes of identical metals (Bp.i = Ep ,2 = Ep) by similar chemical linkage 
(nitrogen-gold affinity in the specific situations discussed below). In this case 

Ai ~ A 2 , (3) 

and Eq. © yields 

£ 0,2 — £ 0,1 ~ EAi — EA 2 + ^im ,2 — ^im,l ■ (4) 

-V-" '-V-' 

from exp. from theory 

Here, the under braces indicate the method to be used below for evaluating the 
corresponding quantities; the LHS can be estimated from available experimental data, 
the RHS can be obtained theoretically via quantum chemical calculations. 

To validate the disentangling scheme proposed here, on which the basic Eq. o 
relies, we will consider molecular junctions based on BDCN and 44BPY. 

4 . 1 . Quantities estimated from experimental data 

The LUMO energy offset for the BDCN molecule can be deduced from the 
experimental value of the transition voltage Vt —>■ Vt,i = 1.69 ± 0.05 V extracted from 
the current-voltage {I — V) curve measured at zero gate potential {Vq = 0); see the 
supplementary information of ref. [2] . Let us briefly remind that the transition voltage 
Vt represents the bias at the minimum of the Eowler-Nordheim quantity ln(//U^) [8]. 
Because the experimental I — V curve [5] turned out to be practically symmetric 
[I[V) ~ —/(—U)], the LUMO energy offset can be estimated as [H] 

£0 = ^Vt, (5) 

which yields [55] 

£ 0,1 = 1.46 ±0.04eV. (6) 

On the other side, the LUMO energy offset for 44BPY-based junctions deduced 
via thermopower data 0121] is 

£ 0,2 = 1-53 ±0.08 eV. (7) 


4-2. Quantities estimated via quantum chemical calculations 

4-2.1. Electron attachment energies The results of the quantum chemical calculations 
for the lowest electron attachment energies EA 12 entering Eq. (0) are collected in 
Table [T] In addition to the values obtained within the EOM-CCSD, EOM-CC2, 
and ADC(2) methods described in Sec. 0 values obtained via energy difference (A- 
) CCSD and DFT methods [27] are also presented there. In the latter, the lowest 
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Method 

EAbdcn (eV) 

EA 44 BBV (eV) 

EAbdcn - EA44Bpy (eV) 

EOM-CCSD 

0.717 

0.032 

0.685 

EOM-CC2 

1.047 

0.360 

0.687 

ADC(2) 

1.107 

0.370 

0.737 

A-CCSD 

0.678 

0.0043 

0.678 

A-DFT 

1.127 

0.444 

0.683 


Table 1. Results for the electron attachment energies EA of the isolated 
molecules (BDCN and 44BPY) considered in this study computed by means of 
various quantum chemical calculations indicated in the left column and described 
in the main text. The geometries of the neutral molecules have been optimized 
at DFT/B3LYP/aug-cc-pVDZ level. Notice that although the absolute values 
EAbdcn and EA44BPY niay significantly depend on the method utilized, their 
difference (last column) is quite insensitive to it. 


electron attachment energy is estimated as the difference between the ground state 
energies {£) of the neutral and anionic species at the equilibrium geometry of the 
neutral molecule (M=CCSD, DFT) 

~ £M.neutral £m. anion- (S) 

The inspection the values given in Table [T] is instructive. It reveals that although the 
absolute values of the electroafhnities EAi = EAbdcn and EA2 = EA44BPY for 
the two molecules considered significantly depend on the quantum chemical method 
utilized, the differences EAi — EA 2 deduced by using the aforementioned methods is 
within the experimental accuracy (c/. Eqs. (jH]) and ([7])). 


4-2.2. Image charge effects An extensive analysis of the effect of image charges in a 
two terminal setup was presented recently [22) . Therefore, only a few details will be 
given here. The interaction energy between two infinite planar electrodes and a point 
charge e located at z in vacuo can be expressed as 


{z) = — i - 2 ' 0 ( 1 ) + '0 


4d 


+ 


Z — Zf 


I _ g-p(2-Zs) 


Zt- z 


I _ g-ni^t-z) 


(9) 


where tpiz) = d lnr(z)/c?z is the digamma function. Eq. ([H]) is obtained from the 
expression deduced within classical electrostatics [56] by inserting the expressions in 
the square parentheses, which ensure that the limits lim^-^^^ ^ remain finite 

and provide good fits of the microscopically calculated potential for the single-plane 
problem {z^Zg, z^zt) [57) [SS]- The positions Zg^t {zs < Zt) of the image planes are 
outwardly shifted by zq from the electrode surfaces z' ^ = Zs,t T zq [Ml HO]: where zq 
represents a quantum correction to the classical result. Numerical values appropriate 
for gold electrodes [Au(lll) faces] are gt ~ 1.25bohr“^ and zq ~ 1.58 A [22] . 

Eq. dni cannot be directly applied to a real molecular junction. Contrary to the 
usual claim [smsiig, for cases relevant for molecular electronics [HITT], realistic 
LUMO’s are not point-like but rather extended over the entire molecule. LUMO 
spatial distributions molecules considered in this study are shown 

in Figures [2] and m Because spatial densities of Kohn-Sham LUMO’s are completely 
unphysical and Hartree-Fock LUMO’s may represent a too crude approximation, like 
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in refs. [Ml HZ], we have calculated the natural orbital expansion of the corresponding 
reduced density matrices at the EOM-CCSD level. For the extra electron, we obtained 
that a single natural orbital almost entirely exhausts the natural orbital expansion; for 
BDCN and 44BPY, the weights are 98.1% and 97.7%, respectively. So, this method is 
indeed best suited to describe the spatial distribution of the extra electron in molecules 
with n-type (LUMO-mediated) conduction. 



Figure 2. The almost singly occupied natural orbital corresponding to the 
anion’s extra electron of the BDCN*“ anion (“LUMO”) obtained via EA-EOM- 
CCSD/aug-cc-pVDZ calculations is delocalized over the whole molecule. The 
LUMO density presented here and in Figure 0 was generated by using Gabedit 

ED- 

As visible in Figures [D and EJ rather than being strongly peaked close to the center 
(which would have justified to assume a point-like LUMO), the natural orbital densities 
pLUMO^y.^ of the extra electron is found to be spread over the whole molecules. 
Therefore, the LUMO energy shift driven by image charges should be calculated by 
appropriately weighting Eq. o 

= r ( 10 ) 

J Zs 

where Pi'^^^(-z) = / dx dy p^^^^ {r) is the LUMO density along the molecular axis 

z. 

For properly estimate the image-driven shifts ^im,i = ^im.BDCN and ^im ,2 = 
^im,44BPY via Eq. ([Toll , attention should be paid to the difference between the 
experimental setups employed in refs. [2] and [6], respectively. This difference is 
illustrated by the two cartoons in the left and right panels of Figure 01 

In an asymmetric STM-setup like that used in the measurements for 44BPY-based 
junctions considered here [B], the usual assumption El 0510] of an infinite plate in the 
immediate vicinity of the (nitrogen) atom at one molecular end is justified only for one 
electrode (STM substrate). As shown recently 05], to model an atomically sharp {e.g., 
pyramidal) STM-tip with a height ndyff , one can consider an infinite planar electrode 
displaced from the tip appex by an effective number of n Au(lll) layers. That is (for 
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Figure 3. The almost singly occupied natural orbital corresponding to the 
anion’s extra electron of the 44BPY*“ anion (“LUMO”) obtained via EA-EOM- 
CCSD/aug-cc-pVDZ calculations is delocalized over the whole molecule. 



Figure 4. Cartoons illustrating the basic features of an asymmetric STM setup 
with a sharp STM-tip [6] (left panel) and the less sharp electrodes of a rather 
symmetric electromigrated junction [2] (right panel). 


atom notation see FigureSl Zg = zn^ — Aau-n + Zq and Zt = Zn^ + (Iau-n — zo + 

{dill — 2.354 A, dAu-N — 2.336) [22]. An estimate 

~-1.39eV (11) 

is obtained by using n Ri 3, a value that turned out to be in excellent agreement with 
the experimental data on 44BPY-based junctions analyzed in ref. [22|. 

In the same spirit, to model the (basically symmetric) experimental setup of the 
electromigrated BDCN-based junctions of ref. [2], we will use Zg = zni —dAu-N + Zo — 
mdiii and Zg = zn 2 + dAu-N — 2o + which amounts to consider image planes 

displaced by an effective number m of Au(lll) layers. Taking a value of m smaller 
than n is in accord to the fact that the two gold electrodes in an electromigrated setup 
|2| are not so sharp as an STM tip |6]. Therefore, to get a simple estimate we will 
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assume m ~ n/2, which yields via Eq. (HOD 

~-0.68eV. (12) 

j.3. Validation of the disentangling scheme 

With the numerical values given by Eqs. ©, ©, (II2D, and m and the first line of 
Table [H the following values of the LHS and RHS of Eq. (g]) are obtained 

£ 0 , 2 -£ 0,1 = 0.07±0.12eV, (13) 

EAi — EA 2 + — $im,l ~ — 0.02 eV. (14) 

So, the values of Eqs. 0 and m are in accord with Eq. o within errors. 
Concerning the estimates of Eqs. 0 and 0, we note that they are not substantially 
affected by the values chosen above for n and m. Since the robustness with respect 
to reasonable changes in n has been analyzed in ref. |22j . we only present here the 
m-dependence (see Figure ED. Based on these results, we conclude that the above 



Figure 5. Dependence on the effective number (m) of extra gold layers of 
the image-driven LUMO energy shift ^im,m in an electromigrated junction, as 
schematically presented in the right panel of Figure [J] See the main text for 
details. 


estimates for and $im ,2 are accurate within ^ 0.1 eV, which is consistent to 

experimental inaccuracies expressed in Eq. 0. 

To end this section, we note that, like for 44BPY-based junctions (see Figure 
4(b) of ref. [l^), gold atoms linked to a BDCN molecule do not substantially affect 
the LUMO spatial distribution. This aspect, which is visible in Figure O is relevant: 
it demonstrates that, for the molecules under consideration, corrections due to image 
charges are not dramatically affected by the cutoff procedure close to electrodes. 


5. Conclusion 

The energetic alignment £q relative to electrodes’ Fermi energy of the dominant 
orbital represents a key quantity that controls the charge transport by tunneling in 
molecular junctions. Disentangling £o in contributions with clear physical origin may 
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be important not only for fundamental nanoscience but also for designing electronic 
nanodevices. 

Methodologically, to validate such disentangling schemes, it is preferable to 
compare two molecular junctions that basically differ in a single respect. The 
disentangling scheme analyzed recently [22] , where we have considered junctions based 
on the same molecule (44BPY) but placed in different environments (solvent [42] versus 
ambient conditions i), was a first step in this direction. As a further step in the 
same vein, in this paper we have considered a disentangling scheme for nanojunctions, 
wherein the two molecules considered (BDCN and 44BPY) are different but their 
chemical linkage to electrodes is similar (nitrogen-gold affinity). 

We believe that the validation of the presently proposed disentangling scheme 
against available experimental data m i is noteworthy. Still, considering more 
transport data for molecular junctions exhibiting n-type conduction to generalize the 
proposed method beyond the two aforementioned cases is highly desirable. I — V 
curves for extended bias ranges well beyond the Ohmic regime (|P|.^V^), allowing to 
determine the transition voltage Vt and thence the MO energy offset Eq [cf Eq. ([5|)], 
supplemented by thermopower data |S| or employing electrodes with different work 
functions [3] as evidence for a (LUMO-mediated, n-)type of conduction, would be best 
suited for this purpose. Oligophenylenes with isocyanide linkages NC-(C 6 H 4 )„-CN, 
i.e. series with several (n) phenylene rings instead of the single (n = 1) ring of the 
BDCN=NC-C 6 H 4 -CN considered above and in ref. |2|, may represent good candidates 
for such investigations. Reliable quantum chemical methods like those used here or 
elsewhere m\ can still be applied for molecular species with up to n = 3 — 4 rings. 
Unfortunately, we were unable to find such experimental transport data in existing 
publications, which are very often restricted to the Ohmic conductance. Still, we 
hope that the present theoretical study will encourage accompanying experimental 
(and further theoretical) efforts to validate similar disentangling scheme that could 
certainly contribute to a better microscopical understanding of the nanotransport. 

We end with the following technical remark. To validate the disentangling scheme 
proposed in ref. [52] we have resorted to A-DFT calculations. As compared to 
more elaborate quantum chemical methods, the A-DFT method is computationally 
considerably less demanding. As revealed by the comparison between the first and 
the last line in Table |T] and also discussed elsewhere [27] , A-DFT-based estimates 
for Eq of a given molecular junction may not be satisfactory. Still, differences 
£ 0,1 ~ £ 0,2 between relevant MO offsets Eo,i and Eo ,2 characterizing different (but 
not too different) molecular junctions estimated within A-DFT can be of an accuracy 
comparable to those based on the computationally very costly EOM-CCSD, which 
represents the state-of-the-art of quantum chemistry of medium size molecules. This 
is also an important aspect, as it allows to understand differences between properties 
of various nano junctions by resorting to lower cost computational approaches. 
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